Ground state of the random-bond spin-1 Heisenberg chain 
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Stochastic series expansion quantum Monte Carlo is used to study the ground state of the an- 
tiferromagnetic spin-1 Heisenberg chain with bond disorder. Typical spin- and string-correlations 
functions behave in accordance with real-space renormalization group predictions for the random- 
singlet phase. The average string-correlation function decays algebraically with an exponent of 
-0.378(6), in very good agreement with the prediction of —(3 — y/E)/2 ~ —0.382, while the average 
spin-correlation function is found to decay with an exponent of about -1, quite different from the 
expected value of -2. By implementing the concept of directed loops for the spin-1 chain we show 
that autocorrelation times can be reduced by up to two orders of magnitude. 

PACS numbers: 75.10.Jm, 75.40.Mg, 75.50.Ee 



I. INTRODUCTION 

The ground state of quantum spin chains displays a 
rich variety of physical phenomena such as quasi-long- 
range order, quantum phase transitions and hidden order 
parameters. Some models are exactly solvable by Bethe 
ansatz techniques, while others have to be approached 
using approximate techniques such as spin-wave theory, 
renormalization groupu(RG) methods, bosonization or 
numerical simulations .El Some methods, such^as the den- 
sity matrix renormalization group technique,!] have been 
developed in the context of spin chains. 

The introduction of disorder into quantum spin chains 
appears to create a very complex problem, yet there 
has been remarkable success in understanding and solv- 
ing certain strongly disordered spin chains. Using a 
real-space renormalization group method Ma, Dasgupta 
and Hu were able to obtain many results for the bond- 
disordered spin- 1/2 Heisenberg chain. a Fisher manage 
to solve the RG equations exactly for the spin- 1/2 case,l 
and he showed that the ground state is the so-called 
random-singlet phase, where all spins pair up and form 
singlets over arbitrarily large distances. The character- 
istic time scale is an exponential function of the length 
scale, resulting in a dynamic critical exponent z = oo. 
The random-singlet phase is gapless and average correla- 
tion functions decay algebraically, while typical correla- 
tions decay exponentially. Due to the singlet coupling be- 
tween spins the decay exponents are predicted to be equal 
for all components of the correlation function, even if the 
Hamiltonian is anisotropic. In the RG calculation any 
initial randomness flows to the infinite-disorder random- 
singlet fixed point. Many of these striku^-prediction 
have been confirmed by numerical studies. BtjQtJ 

The physics of the clean spin-1 chain is dramatically 
different from the spin- 1/2 chain due to the presence of 
the Haldane gapp accompanied by exponentially decreas- 
ing spin-correlation functions and a non-zero string-order 
parameter. The effects of adding disorder to the spin- 
1 chain are more controversial, since the RG equations 
cannot be solved analytically in this case. By mapping 
the spin-1 chain to an effective spin- 1/2 system Hyman 



and Yang argue, that the spin-1 chain is stable against 
weak disorder,E2l but for a sufficiently strong disorder 
the chain undergoes a second order phase transition to 
the random-singlet phase. A density matrix renormaliza=j 
tion group study finds no evidence of such a transition,!!! 
but Hyman and Yang suggest that the disorder was not 
strong enough in the numerical study£3 However, for 
the same disorder distribution a recent quantum Monte 
Carlo (QMC) study fipds evidence of a transition to the 
random-singlet ntLase,li3 and a new RG study also sup- 
ports this result. L3 

In a further attempt to resolve the nature of the ground 
state of the strongly disordered spin-1 chain our goal 
is to study the average and typical spin- and string- 
correlation functions. The decay and distribution of these 
functions is one of the hallmarks of the random-singlet 
phase, and has been studied numerically to cnnfirm the 
random-singlet picture for the spin- 1/2 chainMrEl To ob- 
tain this goal we perform QMC simulations at a tem- 
perature low enough so that the observables in question 
have obtained their ground state expectation values.— We 
use the stochastic series expansion QMC algorithmic in 
our calculations, and by applying the recently introduced 
concept of directed loopaia to the spin-1 chain we demon- 
strate that autocorrelation times can be decreased by up 
to two orders of magnitude in the cases we have exam- 
ined. 

The outline of the paper is as follows: in Sec. II we 
review the basics of the stochastic series expansion. In 
Sec. IH we implement the method of directed loops for 
the spin-1 case, and show the effects on the autocor- 
relation times. The method used to determine average 
and typical correlations is discussed in Sec. IV. Results 



for the ground state behavior of the random-bond spin-1 
chain are shown in Se c. |v[ We conclude with summary 
and discussion in Sec. VI. 



II. OPERATOR-LOOP ALGORITHM 

The stochastic series expansion is described in detail 
elsewhere. I!3 In this section we give a brief summary in 
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order to explain the new features we introduce in Sec. III. 

The model treated in this article is the one-dimensional 
antiferromagnet spin-1 Heisenberg chain with bond dis- 
order. The Hamiltonian, H , is given by, 



H — JiSi ■ Sj+i, 



(1) 



where S denotes a spin-1 operator and the positive cou- 
pling parameter J^, is randomly distributed. 
The partition function, Z, is Taylor expanded, 



(a\H m \a), 



(2) 



a m— 



in the basis {\a)} = {\Sf, S%, .., Sff)}, where N is the 
number of spins. The inverse temperature is denoted by 
P. 

Next we rewrite the Hamiltonian as a sum over diago- 
nal and off-diagonal operators, 



/v 



H 



Jb(H 1 j, 



H- 



2,b) 



(3) 



6=1 



where N is the number of spins and b denotes a bond 
corresponding to a pair of interacting spins j(b) and k(b). 
The operators are given by 



#i,6 — C - Si (1A S, 



(4) 



1 



#2,6 - -{S m S m +S j{b] S k{b) ), 



where C is a constant inserted to assure that the weight 
function is positive for all configurations. 

To simplify the Monte Carlo update we introduce an 
additional unit operator -ffo.o = 1- Inserting the Hamil- 
tonian given by Eq. (^) into Eq. (Q), and truncating the 
sum at m = L we obtain 

z = ^ __!^__ (a | fl J bi H aubi \a), (5) 

a S L ' i=l 

where Sl denotes a sequence of operator-indices 



S L = (oi, 61)1, (a 2 , 6 2 ) 2 , (Oi, &z,)j 



(6) 



with aj = 1,2 and 6; = 1, . . . , AT or (a,, &,) = (0, 0) and n 
is the number of non-unit operators in Si,. 

The Monte Carlo procedure must hence sample the 
space of all states \a), and all sequences Sl- The simula- 
tion starts with some random state \a) and an operator 
string containing only unit operators. One Monte Carlo 
step consists of a diagonal and an off-diagonal update. 
In the diagonal update attempts are made to exchange 
unit and diagonal operators sequentially at each position 
in the operator string. Defining the propagated state 



<p)>~n# ai Ai 



the probability for inserting or deleting a diagonal oper- 
ator at place p in the operator string is given by detailed 
balance, ta 

P([0, 0] p -+ [1, b]p ) = WK°<P-Wy\°(p-1)) (8) 



P([l,6] p ->[0,0] p ) = 



L — n 
L-n+1 



J b Np{a(p-l)\H hb \a(p-l)Y 



The inclusion of J b in the expressions above represents 
the place where the bond disorder appears explicitly 
in the algorithm. Removing or inserting an operator 
changes the expansion power n by ±1. 

The off-diagonal update, also called loop update, is 
carried out with n fixed. Each bond operator H bi = 
Hi tbi + i?2,6i acts only on two spins, Sj/ b A and Sy b A- 
We can therefore rewrite the matrix elements in Eq. (|5|) 
as a product of n terms, called vertices, 



M(a,S L ) = l[W, 



v(i) 1 



(9) 



where the vertex weight W v a) is defined as 

= 

(Sj {b) (p)Sk b) (p)\H bi \S* {b) (p- 1)S| (6) (p- 1)). (10) 

A vertex thus consists of four spins, called the legs of the 
vertex, and an operator. 

The principles of the loop update are: one of the n 
vertices is chosen at random and one of its four legs is 
randomly selected as the entrance leg. The entrance leg 
is given a random state that differs from its initial state. 
One of its four legs is chosen as the exit leg, and its state 
is also changed. In order to satisfy detailed balance the 
exit leg is chosen with a probability proportional to the 
vertex weight W v ^ after the spins have been flipped. 
Thereafter the vertex list is sequentially searched for the 
next vertex that contains the exit spin. This spin be- 
comes the entrance leg of the next vertex and this pro- 
cedure is continued until the loop reaches the original 
entrance leg, when the loop closes. During one Monte 
Carlo step the loop update is repeated until on average 
all the vertices have been updated. 



III. DIRECTED LOOPS FOR SPIN-1 CHAIN 

The operator-loop formulation described above is ap- 
plicable to a wide variety of models O The efficiency of 
the algorithm does, however, depend on the model in 
question. It is particularly efficient for cases where the 
so-called bounce or backtrack process can be avoided. 
If the entrance and exit legs of a vertex coincide, the 
vertex spin configuration is left unchanged and the loop 
backtracks to the previous vertex. This is undesirable 
since the bounce process does not achieve any change 
in the SSE configuration. For models such as the spin- 
1/2 Heisenberg model and spin-1/2 XX model, where the 
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bounce process can be avoided, the autocorrelation times 
are generally significantly shorter than for cases where 
the bounce has to be included. 

Recently Sandvik and Syljuasen showed that, by in- 
troducing the notion of directed loops, it is possible to 
exclude the bounce process for a spin- 1/2 chain in a 
wide parameter regime.Ea The application to more gen- 
eral models is also discussed in Ref. [l6| as well as in 
Ref. [jj. We will now show in some detail how the bounce 
process can be eliminated also for the spin-1 case. 

The detailed balance principle requires that the prob- 
abilities of changing between two state, s and s' , with 
weight functions W(s) and W(s'), satisfy 



W{s)P(s -> s') = W(s')P(s' -> s). 



(11) 



By considering all possible loop configurations that take 
us from configuration s to s', as well as the "time"- 
reversed paths returning to configuration s, Sandvik and 
Syljuasen showed that the detailed balance criteria is sat- 
isfied if the local vertex condition 

W s P(s, h -> a', h) = W s .P{s', l 2 -» s, h) (12) 

is fulfilled. W s denotes the vertex weight, given by 
Eq. p(f), and P(s, l\ — > s', 1%) is the probability to choose 
exit leg I2, given entrance leg l\. The spin configuration of 
the vertex changes from s to s'. This criteria relates the 
transfer probabilities of two processes where the direction 
of the loop segment is reversed, since the entrance and 
exit legs are interchanged. Considering that the probabil- 
ities for choosing different exit legs should sum to unity 
Eq. (O) immediately leads to 



(13) 



We can use Eq. (|12| ) and Eq. (|13| ) to construct the prob- 
abilities, P(s,h — > s',^), needed to perform the loop 
update. It turns out that it is possible to divide the 
space of all processes into subspaces, where only vertices 
in the same group are related to each other by Eq. (12). 
Hence, one may solve the detailed balance equations for 
each transformation group separately, and in many cases 
it is possible to find a solution with no bounce processes. 

We now show how to apply this method to the spin-1 
chain, and how to solve the resulting equations in a man- 
ner that eliminates the bounce process also in this case. 
For the spin-1 model there are in total 52 transformation 
groups, but many of them can be transformed into each 
other by symmetry arguments. There are five irreducible 
groups with two vertices, two groups with three elements, 
and one group with four vertices. The last group is shown 
in Fig. 1, where the operator of each vertex is shown as 
a horizontal line. The operator acts on the two spins 
above it, resulting in the spin configuration shown below 
it. Each row represent all vertices that can be reached 
from a given entrance leg with a given spin projection. 
Eq. (O) gives us the following set of equations 
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FIG. 1: An example of a transformation group with four ver- 
tices. The solid (empty) diamonds represent a spin projection 
of +1 (-1), and the circles the projection of 0. The lines with 
arrows are the directed loop segments. A row in the figure 
shows all possible updates of the vertex for a given entrance 
leg and spin state. The smaller symbols at the beginning and 
end of the loop segment show the entrace and exit spin states 
after the update. 



W 2 = 0,2 + b + e + f 
W 3 = as + c+e + g 
Wi = a± + d + f + g. 

The weight on the left-hand side refers to the weight of 
the bare vertex, given by Eq. (|l0|), while the weights on 
the right-hand side refer to the product of vertex weight 
and the transfer probability. Processes that are related to 
each other by changing the direction of the loop segment, 
and interchanging spin configurations, are given equal 
weights, according to Eq. (|12|). It is possible to solve 
this set of equations under the restriction that all bounce 
processes cii are zero. There is not a unique solution 
since the number of unknowns exceeds the number of 
equations, and we have chosen the solution which gives 
all the allowed vertices equal weight. 

For the spin-1 case the groups can be classified accord- 
ing to the change in spin projection on the entrance spin. 
All the groups where the spin projection changes by one 
can be solved so that the bounce is eliminated. There 
are, however, also four groups where the spin projection 
changes by two, shown in Fig. ^. The bounce cannot be 
eliminated in this case. The equations for the group in 
the bottom right corner are of the form 



Wi 



at 



(14) 



Wi 



a 1 



C+l 



(15) 
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FIG. 2: The four transformation groups for the spin-1 chain, 
in which the magnetization changes by two. The solid (empty) 
diamonds represent a spin projection of +1 (-1), and the cir- 
cles a projection of 0. 
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and the bounce processes a\ and et2 cannot be eliminated. 
We can, however, solve all groups shown in Fig. || so that 
the only processes allowed are the bounces. This means 
that any attempt to change the magnetization by two 
will lead to an immediate bounce and the loop will close. 
In practice this idea can be implemented by attempting 
to change the magnetization only by one, and accepting 
50% of the attempts when the initial spin state is ±1. 

To study the efficiency of the new loop move, inte- 
grated autocorrelation times have been calculated. The 
normalized autocorrelation function is defined as, 



A Q {t) = 



(Q(i + t)Q(i))-(Q(i))< 

m?) - <q«> 2 



where i and t are measured in units of Monte Carlo steps, 
with one Monte Carlo step defined as one diagonal up- 
date and one loop update. The integrated autocorrela- 
tion time is defined according to 



Tint [ 



1 OO 



(17) 



4=1 



In Fig. m integrated autocorrelation times are shown for 
calculations using the old and the new loop updating 
method. The calculations are done on a system con- 
sisting of 32 spins without any disorder in the coupling 
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FIG. 3: Integrated autocorrelation times as a function of in- 
verse temperature, /3J, for a system with 32 spins without 
bond disorder. The open symbols represents the directed loop 
updating method, where bounces have been excluded and the 
solid symbols are calculated using the original loop move. 



parameters. The autocorrelation times are calculated for 
the energy and for the spin- and string-correlation func- 
tion, which we now define. 

The spin-correlation function C Z (N) for a chain with 
N spins is measured at the greatest distance between 
two spins, which for periodic boundary conditions equals 
N/2, 



C Z {N) 



1 

N 



N 



(18) 



i=i 



and the string-correlation function O z (N) is similarly de- 
fined by, 



1 N 

° Z ( N ) = n ^ {S: ex P^ 5 + ^+2 



i=l 



+s z 



}S 



(16) The energy is calculated according tct_ 



E 



1 



(19) 



(20) 



where n is the number of non-unit operators in the op- 
erator string. The energy autocorrelation time is least 
affected by the new loop update, as is clearly seen in 
Fig. H The reason is that the energy is not directly de- 
pendent on the loop move, since the energy is given by 
the number of operators which is determined in the di- 
agonal update. The exclusion of the bounce shortens the 
autocorrelation time by up to two orders of magnitude 
for the spin correlation at low temperatures. 
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IV. DETERMINING AVERAGE AND TYPICAL 
CORRELATIONS 

To calculate ground state expectation values in 
strongly disordered systems using QMC is not an easy 
task. The standard deviation in the disorder averaged 
expectation values has two sources: the standard statis- 
tical fluctuations, and the sample to sample variation. In 
strongly disordered systems the latter is often dominant 
and it is therefore advantageous to perform short simula- 
tions of a large number of disorder configurations. This is 
problematic for two reasons. First, each sample requires 
an adequate equilibration time, which becomes a very 
significant part of the simulation, particularly for short 
simulations. Second, it may be hard to determine at 
what temperature the expectation values have converged 
to their ground state values. Repeating the calculations 
at many temperatures is very time demanding. 

Recently^ technique was proposed to alleviate these 
difficulties .El The basic idea is to perform a very small 
number of Monte Carlo steps at a series of decreasing 
temperatures. The method builds on two realizations. 
First, given an equilibrium Monte Carlo configuration at 
inverse temperature /3 it is possible to construct an ap- 
proximate equilibrium configuration at inverse tempera- 
ture 2(3. Second, due to the extremely short autocorre- 
lation times equilibrium is reached quickly. Due to the 
small number of Monte Carlo steps necessary at each 
temperature the methods allows very rapid calculation 
of disorder averages down to very low temperatures. 

Now we will briefly describe the implementation of the 
method. The energy expectation value, see Eq. (pOf), is 
proportional to the length of the operator string and the 
temperature. This means that the length of the opera- 
tor string is proportional to the inverse temperature (3 as 
the ground state is approached. The calculation begins 
at a high temperature, where a few equilibration Monte 
Carlo steps are performed, followed by a few steps when 
measurements are carried out. Thereafter the tempera- 
ture is halved at the same time as the operator string 
is doubled. The procedure is repeated and results in a 
series of measurements at decreasing temperatures. This 
is carried on until disorder averages of the expectation 
values show no further temperature dependence. 

In the random-singlet phase the spin-correlation func- 
tions have a very broad distribution. The average of the 
correlation functions will be dominated by very strong 
correlations, while the typical correlations are much 
weaker. In order to characterize this behavior we measure 
not only the average of the correlation function, but also 
the average of the logarithm of the correlation function, 
which defines the typical correlation function, 

1 N 

\ogC? yp (N) = - 5>g<S?S^ f ). (21) 

i=l 

The method described above is suitable for determin- 
ing average correlation functions. However, it cannot 



be used to accurately determine typical correlation func- 
tions. The reason is that small values of the correlation 
functions give large contributions to the typical correla- 
tions. Therefore, the whole distribution of the correla- 
tion function needs to be determined quite accurately to 
investigate the typical behavior. In the above method 
the statistical errors for a single disorder configuration 
are very large, due to the small number of Monte Carlo 
steps. Hence the weak correlations, which are impor- 
tant for the typical correlations, drown in the statistical 
noise. In particular, many negative values will have to 
be discarded since the logarithm is defined only for posi- 
tive arguments. In order to determine typical correlation 
functions we have therefore done the simulation in two 
stages. First we use the above method to determine the 
temperature at which average correlation functions have 
converged in temperature. Thereafter we perform much 
longer simulations at this temperature to determine the 
typical correlation functions. 

V. RESULTS FOR THE RANDOM-BOND 
CHAIN 

The properties of the random-singlet phase are pre- 
dicted to be independent of the underlying Hamiltonian, 
as long as the ground state of two nearest- neighbor spins 
is a singlet. In particular, the decay exponents of the spin 
and string-correlation functions-should be the same for 
the spin-1 and spin-1/2 chain.td In the random-singlet 
phase it is predicted that the average correlation func- 
tions decay algebraically,[3 

C z (N)otN- 2 (22) 

for large AT, and 

O z {N) ocN-^ ~N-°- 382 , (23) 

while typical correlations decay with a stretched expo- 
nential as 

C* typ {N) cx exp(-AViV), (24) 

where A is a non-universal constant. 

To make sure that the QMC method is capable of cor- 
rectly determining the correlation functions for a strongly 
disordered spin system we have first performed a calcula- 
tion on the spin-1/2 XX chain with bond disorder. The 
XX chain is described by the following Hamiltonian, 

H = J2Ji(sfsf +1 + sysy +1 ), (25) 

i 

where S denotes a spin-1/2 operator. The ground state 
expectation values of the correlation functions for this 
system have been calculated by exact diagonalization us-, 
ing the Jordan- Wigner transformation to free fermions.El 
We compare QMC results for the average correlation 
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FIG. 4: A comparison between exact diagonalization and 
QMC results for the correlation functions of the spin- 1/2 XX 
chain with bond disorder. Lines show algebraic decay with ex- 
ponents predicted from RG calculations. The lines are drawn 
through the last exact diagonalization data point. 



FIG. 6: The spin- and string-correlation functions for the 
bond-disordered spin-1 chain. Disorder averages are taken 
over at least 1000 configurations. The solid lines are linear 
fits to the data points. 
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FIG. 5: The temperature dependence of the string- and spin- 
correlation function for a system with 64 spins. 



functions with the diagonalization data in Fig. || and 
there is very good agreement. The lines are drawn 
through the last diagonalization data point with the pre- 
dicted slopes of -2 for the spin-correlation function and 
-0.382 for the string-correlation function. We can see 
that already for these relatively small system sizes the 
data displays almost linear behavior with slopes close to 
the predicted values. The finite-size corrections are more 
pronounced for the string-correlation function. 

Next we consider the spin-1 chain. Recent numerical 
studies of the disordered spin-1 chain have used bonds 



distributed according to a box distribution 

P(J) = { l/W fOT 1 _ ^ " Ji ~ 1 + 
I otherwise. 



w 

2 



(26) 



A density matrix renormalization group studji did not 
find a transition to the random-singlet phase,til even fps 
the strongest disorder (W = 2), while a QMC studyfc 3 ! 
found evidence for a transition around W — 1.8 p-, 1.9. 
According to the RG analysis by Hyman and Yangta the 
phase transition occurs for the power law distribution, 
P(J) oc J~ 33 , which represents stronger disorcLei than 
the box distribution. Finally, a recent RG studyEj finds 
that the transition occurs exactly at W = 2. Close to 
the critical point the RG studies find that the Haldane 
phase is Griffiths-like, with no gap, finite string order and 
exponentially decaying spin-correlation functions. 

In this work we do not attempt to find the critical 
bond distribution. Instead we try to resolve the con- 
troversy concerning the box distribution. We focus on 
the widest possible box distribution and make an exten- 
sive QMC simulation of the correlation functions to de- 
termine whether they behave according to the random- 
singlet predictions. 

The temperature dependence of the average string- and 
spin-correlation function for system size N = 64 is shown 
in Fig. ^. At each temperature we performed 20 equili- 
bration steps followed by 60 measurement steps. The 
results seem to have converged (within statistical error 
bars) at approximately /3 w 32000. The convergence 
temperatures for different system sizes are presented in 
Table |. At least 1000 disorder distributions are used in 
the calculations. For these system sizes the convergence 
temperature is given by (3 = iV 3 /8. 

In the second stage of the simulations 10 3 calibration 
steps are done, followed by at least 2 x 10 3 measuring 
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FIG. 7: The typical spin- and string-correlation functions as 
a function of the square root of the distance. The solid lines 
are linear fits to the data points. 



steps. Averages are calculated over 1000 disorder con- 
figurations. Results for average and typical correlation 
functions are shown in Fig. ^ and Fig. fj] respectively. 

Straight lines are adjusted to the data points to deter- 
mine the decay exponents in Fig. ^. In the log-log plot 
the string-correlation function displays very linear behav- 
ior, indicating that the function does decay algebraically. 
The decay exponent is found to be -0.378(6), in excel- 
lent agreement with the random-singlet prediction of - 
0.382. The spin-correlation function does not adjust to a 
straight line as well as the string correlation. The gradi- 
ent of the curve increases with increasing system size. If 
a line is adjusted it results in a decay exponent of about 
-1.04(1), quite different from the predicted value of -2. 
The curvature of the function does indicates that we can- 
not see the true decay exponent with the limited system 
sizes used, and we cannot rule out the possibility that the 
exponent will converge to -2 for much larger system sizes. 
It is, however, remarkable that the spin-correlation func- 
tion should display such dramatic system size dependence 
while the string-correlation functions has converged. 

Next we consider the behavior of the typical correlation 
functions in Fig. 0. Straight lines should result if we 
plot the logarithm of the typical correlation functions 
versus the square root of the systems size. The QMG 



System size 


Inverse temperatures /3J 


8 


64 


16 


512 


32 


4096 


64 


32768 



TABLE I: Inverse temperatures needed to reach ground state 
properties for different system sizes. 



FIG. 8: Distribution of the spin-correlation function. 
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FIG. 9: The scaled distribution function for the spin- 
correlation function. 



data lies close to a straight line, but there does seem 
to be a slight systematic curvature upwards. This can 
be explained by the difficulties in measuring the typical 
functions. A small percentage of the measurements are 
discarded since they are negative and this will result in 
too large an estimate. Since more points are dropped for 
the larger system sizes this should result in an upward 
curvature. For the 32-site system about 0.1 % of the 
points were dropped, while about 0.6 % of the points 
were dropped for the 64-site system. We believe that this 
explains the slight systematic deviation from the straight 
line. 

We also considered the distribution of the correlation 
functions. According to the RG predictions the distribu- 
tion of 



N~ 



(27) 



should scale to a fixed distribution with /i = 2. In Fig. 
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FIG. 10: Distribution of the string-correlation function. 



we show the distribution of the logarithm of the spin- 
correlation function. The distribution becomes broader 
and the peak shifts to the left with increasing N. The 
best scaling plot was obtained with fi — 2.2, shown in 
Fig. H The agreement with the RG predictions is again 
very good, but we were not able to scale the distribution 
of the string-correlation function. In Fig. |l^ the unsealed 
distribution is shown, and as can be seen the distribution 
does broaden with increasing N, but much less than for 
the spin correlation function. In order to scale the distri- 
bution so that the peaks line up, a large value of fx « 10 
is needed, but then the rest of the distribution does not 
scale well. 

Finally, we also studied the distribution of the local 
susceptibility, 



Xloc — 



1 N 

-Y 



dr(S!(0)S!(T)), 



(28) 



which is also predicted to scale according to Eq. ( |27[ ). 
In agreement with the earlier QMC study performed at 
higher temperatures we also obtain the best scaling for 
fx = 2.5. 



VI. SUMMARY AND DISCUSSION 

We have done extensive QMC simulations of the an- 
tiferromagnetic spin-1 Heisenberg chain with a uniform 
bond distribution extending all the way to zero bond 
strength. The calculations were performed at a tempera- 
ture low enough for all observables to obtain their ground 
state expectation values. The decay exponent for the 
string-correlation function, as well as the typical correla- 
tion functions, are in good agreement with renormaliza- 
tion group predictions for the random-singlet phase. Our 
results for the average spin correlation function yields 
an exponent of -1.04(1) which differs significantly from 



the RG prediction of -2. However, this result has not 
converged, but does display some finite-size effects. The 
distribution of the spin-correlation function is well de- 
scribed by RG predictions, but the distribution of the 
string-correlation function did not scale as predicted. 

The very clear algebraic decay of the string order 
strongly indicates that the system in no longer in the 
Haldane phase. The exponent we obtained (-0.378(6)) 
is in excellent agreement with the prediction for the 
random-singlet— phase (-0.382), and very different from 
the predictionE2l for the critical point separating the 
random singlet phase from the Haldane phase (-2/3). 
We therefore believe that our results are in agreement 
with earliet, Monte Carlo resultstj and a recent RG 
calculationE-3 showing that the box distribution extend- 
ing all the way to zero is enough to drive the system to an 
infinite-randomness fixed point. The deviations from the 
random-singlet predictions could be explained by strong 
finite-size corrections. An alternative scenario is that the 
strong randomness fixed point for the spin-1 chain has 
slightly different properties from the spin- 1/2 random- 
singlet fixed point. 

To numerically determine the properties of the 
random-singlet phase for the strongly disordered spin-1 
chain is much harder than for the spin- 1/2 chain since a 
strong, finite disorder is needed to drive the chain to the 
random-singlet phase, according to the RG predictions. 
Disagreement between earlier numerical studies indicate 
the difficulties. The widest box distribution is expected 
to be quite close to the transition to the Haldane state, 
and it is possible that strong finite-size corrections appear 
in certain quantities. It would therefore be of great inter- 
est to study larger systems with even stronger, power-law 
distribution of the bonds. This is unfortunately not an 
easy task to do numerically. Density matrix renormal- 
ization group calculations yield ground state expectation 
values, but the method is most accurate for open bound- 
ary conditions for which finite size effects are stronger. 
Furthermore, it is hard to keep enough states to obtain 
accurate results for large systems with strong bond dis- 
order. Using loop quantum Monte Carlo algorithms it 
is possible to reach fairly large system sizes at quite low 
temperatures, but to reach temperatures low enough to 
accurately determine the decay of correlation functions 
makes the simulations very time consuming. Further- 
more, the expectation value of the correlation function 
becomes very small for large system sizes and long simu- 
lations are needed in order to reduce the statistical errors 
to the required level. 

In conclusion, we believe thaj_our results support 
earlier Mcpte Carlo simulationsEj and a recent RG 
calculationE-3 indicating that the box distribution extend- 
ing all the way to zero is wide enough to reach the 
random-singlet phase. With our present data we are 
not able to determine whether observed deviations from 
predictions for the random-singlet fixed point, most no- 
tably in the decay of the spin-correlation function, are 
due to strong finite-size corrections or slightly different 
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properties of the spin-1 random-singlet fixed point. To 
resolve this issue requires significant computational ef- 
forts and we leave it for future investigations. To obtain 
the present results we have implemented the concept of 
directed loops for the spin-1 chain and demonstrate that 
autocorrelation times are decreased by up to two orders 
of magnitude. 
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